#include "libmesh/quadrature.h"
#include "Assembly.h"
#include "EntropyJumpIndicator.h"
#include "CFDProblem.h"

template<>
InputParameters validParams<EntropyJumpIndicator>()
{
  InputParameters params = validParams<FluxJumpIndicator>();
  return params;
}

EntropyJumpIndicator::EntropyJumpIndicator(const InputParameters &parameters) :
	FluxJumpIndicator(parameters)
{
}

void EntropyJumpIndicator::computeIndicator()
{
	Real sum = 0;
	for (_qp=0; _qp<_qrule->n_points(); _qp++)
		sum += _JxW[_qp]*_coord[_qp]*computeQpIntegral();

	sum /= _current_side_volume;
	Real scale = getParam<Real>("scale");
	sum *= scale;

	{
		Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
		_solution.add(_field_var.nodalDofIndex(), fabs(sum));
		_solution.add(_field_var.nodalDofIndexNeighbor(), fabs(sum));
	}
}

Real EntropyJumpIndicator::computeQpIntegral()
{
	for (size_t i = 0; i < _uh.size(); ++i)
	{
		_cfd_data.uh[i] = (*_uh[i])[_qp];
		_cfd_data_neighbor.uh[i] = (*_uh_neighbor[i])[_qp];
	}
	_cfd_data.reinitInviscous();
	_cfd_data_neighbor.reinitInviscous();
	Real vn = (_cfd_data.vel+_cfd_data_neighbor.vel)*_normals[_qp] + _cfd_data.c+_cfd_data_neighbor.c;

	return (_cfd_data.s - _cfd_data_neighbor.s)/1*vn;
}

void EntropyJumpIndicator::finalize()
{

}

